Integration
What if we wanted to look at population in 20 years given an initial
condition
Two options
explicit solution to differential equation is known; e.g. you can
integrate both sides of the equation! Not always possible but lets look
at a case where it is possible
must be solved by iteration; this is what we do when we can’t
integrate both sides
Explicit Solution is available
source("../R/exppop.R")
exppop
## function (T, P0, r, K)
## {
## P = P0 * exp(r * T)
## if (P > K) {
## P = K
## }
## return(P)
## }
# gives population after any time given an initial population
# 20 rabbits, growth rate of 0.01 how many in 30 years
exppop(T=30, P0=20, r=0.01, K=1000)
## [1] 26.99718
# if we want to see how population evolves over time - generate a time series by running our model for each point in time
initialrabbits = 20
years = seq(from=1, to=100, by=2)
Ptime = years %>% map_dbl(~exppop( P0=initialrabbits, r=0.01, K=1000, T=.x))
# keep track of what times we ran
Ptime = data.frame(P=Ptime, years=years)
ggplot(Ptime, aes(years,P))+geom_point()+labs(x="years",y="Rabbit Population")
Solving using numeric integration
Using a solver….when you can’t do the integration by hand :)
For example, if you added a carrying capacity as threshold where it
stops growing
In this case we integrate by iteration - approximates
analytic integration
Numerical integration with ODE
Implement the differential equation as a function that
- returns the derivative (as a list)
*inputs time, the variable(s) and a parameter list
(it needs time even though you don’t use it)
My convention: name derivative functions starting with d to
remind myself that they are computing a derivative
ODE
Only works for Ordinary Differential Equations - single independent
variable (in our case time)
Partial differential equations - more than 1 independent variable
(e.g x and y if changing in space)
R has a solver called ODE for solving ordinary differential
equations from package desolve
ODE requires
- initial conditions
- values of independent where we want values of dependent variable
(e.g times where we want population)
- the derivative as a function
- a list that contains alll parameter values (or if you have only one
parameter then you can use a single value)
source("../R/dexppop.R")
dexppop
## function (time, P, r)
## {
## dexpop = r * P
## return(list(dexpop))
## }
library(deSolve)
initialrabbits = 20
years = seq(from=1, to=100, by=2)
# run the solver
Ptime = ode(y=initialrabbits, times=years, func=dexppop,parms=c(0.01))
head(Ptime)
## time 1
## [1,] 1 20.00000
## [2,] 3 20.40404
## [3,] 5 20.81623
## [4,] 7 21.23675
## [5,] 9 21.66576
## [6,] 11 22.10344
colnames(Ptime)=c("year","P")
# notice that there are additional pieces of information year, including the method used for integration
attributes(Ptime)
## $dim
## [1] 50 2
##
## $dimnames
## $dimnames[[1]]
## NULL
##
## $dimnames[[2]]
## [1] "year" "P"
##
##
## $istate
## [1] 2 52 105 NA 6 6 0 36 21 NA NA NA NA 0 1 1 NA NA NA
## [20] NA NA
##
## $rstate
## [1] 2.00000 2.00000 99.98839 0.00000 1.00000
##
## $lengthvar
## [1] 1
##
## $class
## [1] "deSolve" "matrix"
##
## $type
## [1] "lsoda"
# this also means you need to extract just the data frame for plotting
ggplot(as.data.frame(Ptime), aes(year,P))+geom_point()+labs(y="Population", "years")
# this also works (of course function can be by order)
Ptime=ode(initialrabbits, years, dexppop,0.03)
colnames(Ptime)=c("year","P")
ggplot(as.data.frame(Ptime), aes(year,P))+geom_point()+labs(y="Population", "years")
additional parameters
You can play a bit with changing your function to something that you
can’t integrate “by hand”
BUT we might want more parameters
to work with ODE, parameters must all be input as a single list;
similar to how we return multiple outputs from a function
see example below..lets add a carrying capacity
R code with carrying capacity
source("../R/dexppop_play.R")
dexppop_play
## function (time, P, parms)
## {
## dexpop = parms$r * P
## dexpop = ifelse(P > parms$carry_capacity, 0, dexpop)
## return(list(dexpop))
## }
# create parameter list
initalrabbits=2
newparms = list(r=0.03, carry_capacity=300)
#apply solver
results=ode(initialrabbits, years, dexppop_play,newparms)
head(results)
## time 1
## [1,] 1 20.00000
## [2,] 3 21.23675
## [3,] 5 22.54993
## [4,] 7 23.94435
## [5,] 9 25.42499
## [6,] 11 26.99719
# add more meaningful names
colnames(results)=c("year","P")
# plot
ggplot(as.data.frame(results), aes(year,P))+geom_point()+labs(y="Population", "years")
# try again with different parameters
alternativeparms = list(r=0.04, carry_capacity=500)
results2=ode(initialrabbits, years, dexppop_play,alternativeparms)
## DLSODA- At current T (=R1), MXSTEP (=I1) steps
## taken on this call before reaching TOUT
## In above message, I1 = 5000
##
## In above message, R1 = 81.5108
##
## Warning in lsoda(y, times, func, parms, ...): an excessive amount of work (>
## maxsteps ) was done, but integration was not successful - increase maxsteps
## Warning in lsoda(y, times, func, parms, ...): Returning early. Results are
## accurate, as far as they go
# look at results
head(results2)
## time 1
## [1,] 1 20.00000
## [2,] 3 21.66575
## [3,] 5 23.47022
## [4,] 7 25.42499
## [5,] 9 27.54257
## [6,] 11 29.83650
colnames(results2)=c("year","P_parm2")
# plot
ggplot(as.data.frame(results2), aes(year,P_parm2))+geom_point()+labs(y="Population", "years")
# compare by combining into a single data frame
both = inner_join(as.data.frame(results), as.data.frame(results2))
## Joining, by = "year"
both_p = both %>% gather(key="model", value="P", -year)
ggplot(both_p, aes(year,P, col=model))+geom_point()+labs(y="Population", "years")
# try playing on your own - modify the function in some way
Difference Equations
What is ODE doing? (iterating in ‘smart ways’)
Similar to “difference equations”
Difference Equations
Population models can be discrete (rather than continuous) So we
could implement them as difference equations and iterate
source("../R/discrete_logistic_pop.R")
# notice how a for loop is used to iterate
# how many rabbits after 50 years given a growth of 0.1
# starting with 1 rabbit - but a carrying capcity of 500
discrete_logistic_pop
## function (P0, r, K, T = 10)
## {
## pop = P0
## for (i in 1:T) {
## pop = pop + r * pop
## pop = ifelse(pop > K, K, pop)
## }
## return(pop)
## }
discrete_logistic_pop(P0=1, r=0.05, K=200, T=50)
## [1] 11.4674
# save results
discrete_result = discrete_logistic_pop(P0=1, r=0.05, K=200, T=50)
# lets also keep the parameters for use later
P0=1
r=0.05
K=200
T=50
Differential Equation, Difference (Iteration by hand)
comparison
Remember we have 3 ways now to calculate population
analytical solution - based on integration (exppop.R) BEST
using an ode solver for numerical approximation (exppop_play.R)
numerical integration using in discrete steps
(discrete_logistic_pop.R)
source("../R/exppop.R")
exppop(P0=P0, r=r, K=K, T=T)
## [1] 12.18249
analytic_result = exppop(P0=P0, r=r, K=K, T=T)
analytic_result
## [1] 12.18249
## [1] 11.4674
# why are they different
# look at trajectories
growth_result = data.frame(time=seq(from=1,to=100))
growth_result$Panalytic = growth_result$time %>% map_dbl(~exppop( P0=1,r=0.05, K=200,T=.x ))
growth_result$Pdiscrete = growth_result$time %>% map_dbl(~discrete_logistic_pop( P0=1,r=0.05, K=200,T=.x ))
tmp = growth_result %>% gather(key="Ptype",value="P",-time)
ggplot(tmp, aes(time,P, col=Ptype))+geom_point()
# try running them for longer time periods to see what happens
# change the value of r, K , P0 - see how it effects the results
Compare analytical, difference and ODE
Finally look at continuous derivative using ODE solve Needs
initial condtions differential equation *parameters
source("../R/dexppop_play.R")
dexppop_play
## function (time, P, parms)
## {
## dexpop = parms$r * P
## dexpop = ifelse(P > parms$carry_capacity, 0, dexpop)
## return(list(dexpop))
## }
# set up using the same parameters
pcompare = list(r=r,carry_capacity=K)
# now run our ODE solver
result = ode(P0, growth_result$time, dexppop_play, pcompare)
head(result)
## time 1
## [1,] 1 1.000000
## [2,] 2 1.051273
## [3,] 3 1.105172
## [4,] 4 1.161835
## [5,] 5 1.221404
## [6,] 6 1.284027
# we already have time - so just extract population
growth_result$Pdifferential = result[,2]
# comapre all 3 approaches
tmp = growth_result %>% gather(key="Ptype",value="P",-time)
ggplot(tmp, aes(time,P, col=Ptype))+geom_point()
# notice Pdifferential is closer to Panalytic than Pdiscrete
Other Examples
All differential and difference equations are approximations The
analytical solution is exact
Notice that differential equations is a bit more accurate!
Diffusion in Theory
Diffusion in 1 dimension
Diffusion in one dimension through time
Diffusion data structure
